function [G_final, S_final, F_final, nIter_final, objhistory_final] = DRCC_Multi(X, rowK, colK, rowW, colW, options, G, S, F)
% Dual Regularized Co-Clustering (DRCC) with
%          multiplicative update
%
% Input
%     X: (mFea x nSmp) data matrix
%         mFea  ... number of words (vocabulary size)
%         nSmp  ... number of documents
%     rowK: number of row hidden factors
%     colK: number of column hidden factors
%     rowW: weight matrix of the affinity graph on features
%     colW: weight matrix of the affinity graph on samples
%     options: Structure holding all settings
%         options.rowLamda ... the regularization parameter. [default: 100]
%         options.colLamda ... the regularization parameter. [default: 100]
%
% You only need to provide the above four inputs.
%
% Optimization problem
%     min ||X - G S F'||^2 + rowLamda tr(G' Lg G)  + colLamda tr(F' Lf F)
%     s.t. G, S, F >=0.
%
% References:
% [1] Quanquan Gu, Jie Zhou. "Co-Clustering on Manifolds", KDD, 2009.
%
%
%   Written by Liang Du (csliangdu AT gmail.com)
%
%

differror = options.error;
maxIter = options.maxIter;
nRepeat = options.nRepeat;
minIter = options.minIter - 1;
if ~isempty(maxIter) && maxIter < minIter
    minIter = maxIter;
end
meanFitRatio = options.meanFitRatio;

rowLamda = options.rowLamda;
colLamda = options.colLamda;

[nDim,nSmp]=size(X);

Lg = ConstructLaplacian(rowW, rowLamda, options);
Lf = ConstructLaplacian(colW, colLamda, options);

selectInit = 1;
if isempty(G)
    G = abs(rand(nDim, rowK));
    S = abs(rand(rowK, colK));
    F = abs(rand(nSmp, colK));
else
    nRepeat = 1;
end


if nRepeat == 1
    selectInit = 0;
    minIter = 0;
    if isempty(maxIter)
        objhistory = CalculateObj(X, G, S, F, Lg, Lf);
        meanFit = objhistory*10;
    else
        if isfield(options,'Converge') && options.Converge
            objhistory = CalculateObj(X, G, S, F, Lg, Lf);
        end
    end
else
    if isfield(options,'Converge') && options.Converge
        error('Not implemented!');
    end
end



tryNo = 0;
nIter = 0;
while tryNo < nRepeat
    tryNo = tryNo+1;
    maxErr = 1;
    while(maxErr > differror)
        % *************************************
        % Update S
        % *************************************
        S = inv(G' * G) * (G' * X * F) * inv(F' * F);
        
        % *************************************
        % Update G
        % *************************************
        tmp = F * S';
        P = X * tmp;
        Q = tmp' * tmp;
        t1 = posX(P) + G * negX(Q);
        t2 = negX(P) + G * posX(Q);
        if rowLamda > 0
            t1 = t1 + negX(Lg) * G;
            t2 = t2 + posX(Lg) * G;
        end
        G = G .* sqrt(t1./ max(t2, eps));
        G = (t1 > 0) .* G;
        clear P Q tmp t1 t2;
        
        % *************************************
        % Update F
        % *************************************
        tmp = G * S;
        A = X' * tmp;
        B = tmp' * tmp;
        t1 =  posX(A) + F * negX(B);
        t2 =  negX(A) + F * posX(B);
        if colLamda > 0
            t1 = t1 + negX(Lf) * F;
            t2 = t2 + posX(Lf) * F;
        end
        F = F .* sqrt(t1./ max(t2, eps));
        F = (t1 > 0) .* F;
        clear A B tmp t1 t2;
        
        
        nIter = nIter + 1;
        if nIter > minIter
            if selectInit
                objhistory = CalculateObj(X, G, S, F, Lg, Lf);
                maxErr = 0;
            else
                if isempty(maxIter)
                    newobj = CalculateObj(X, G, S, F, Lg, Lf);
                    objhistory = [objhistory; newobj]; %#ok<AGROW>
                    meanFit = meanFitRatio*meanFit + (1-meanFitRatio)*newobj;
                    maxErr = (meanFit-newobj)/meanFit;
                else
                    if isfield(options,'Converge') && options.Converge
                        newobj = CalculateObj(X, G, S, F, Lg, Lf);
                        objhistory = [objhistory; newobj]; %#ok<AGROW>
                    end
                    maxErr = 1;
                    if nIter >= maxIter
                        maxErr = 0;
                        if isfield(options,'Converge') && options.Converge
                        else
                            objhistory = 0;
                        end
                    end
                end
            end
        end
    end
    
    if tryNo == 1
        G_final = G;
        S_final = S;
        F_final = F;
        nIter_final = nIter;
        objhistory_final = objhistory;
    else
        if objhistory(end) < objhistory_final(end)
            G_final = G;
            S_final = S;
            F_final = F;
            nIter_final = nIter;
            objhistory_final = objhistory;
        end
    end
    
    if selectInit
        if tryNo < nRepeat
            %re-start
            G = abs(rand(nDim, rowK));
            S = abs(rand(rowK, colK));
            F = abs(rand(nSmp, colK));
            nIter = 0;
        else
            tryNo = tryNo - 1;
            nIter = minIter+1;
            selectInit = 0;
            G = G_final;
            S = S_final;
            F = F_final;
            objhistory = objhistory_final;
            meanFit = objhistory*10;
        end
    end
end

end

function [obj, dV] = CalculateObj(X, G, S, F, Lg, Lf)
MAXARRAY = 500*1024*1024/8; % 500M. You can modify this number based on your machine's computational power.
if ~exist('Lg','var')
    Lg = 0;
end
if ~exist('Lf','var')
    Lf = 0;
end
dV = [];
nSmp = size(X,2);
mn = numel(X);
nBlock = ceil(mn/MAXARRAY);

if mn < MAXARRAY
    dX = X - G * S * F';
    obj_NMF = sum(sum(dX.^2));
else
    obj_NMF = 0;
    Gs = G * S;
    PatchSize = ceil(nSmp/nBlock);
    for i = 1:nBlock
        if i*PatchSize > nSmp
            smpIdx = (i-1)*PatchSize+1:nSmp;
        else
            smpIdx = (i-1)*PatchSize+1:i*PatchSize;
        end
        dX = X(:,smpIdx) -  Gs * F(smpIdx,:)';
        obj_NMF = obj_NMF + sum(sum(dX.^2));
    end
end

if isempty(Lg)
    obj_Lg = 0;
else
    obj_Lg = sum(sum((Lg * G).* G));
end

if isempty(Lf)
    obj_Lf = 0;
else
    obj_Lf = sum(sum((Lf * F).* F));
end
obj = obj_NMF + obj_Lg + obj_Lf;
end

function X = posX(X)
X = (abs(X) + X)/2;
end

function X = negX(X)
X = (abs(X) - X)/2;
end

function L = ConstructLaplacian(W, lamda, options)
nSmp = size(W,1);
if lamda > 0
    W = lamda * W;
    DCol = full(sum(W,2));
    D = spdiags(DCol,0,nSmp,nSmp);
    L = D - W;
    if isfield(options,'NormW') && options.NormW
        D_mhalf = spdiags(DCol.^-.5,0,nSmp,nSmp) ;
        L = D_mhalf*L*D_mhalf;
    end
else
    L = [];
end
end